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Abstract 

It is demonstrated that the minimization of the free energy func- 
tional for hard spheres and hard disks yields the result that excited 
granular materials under gravity segregate not only in the widely 
known "Brazil nut" fashion, i.e. with the larger particles rising to 
the top, but also in reverse "Brazil nut" fashion. Specifically, the lo- 
cal density approximation is used to investigate the crossover between 
the two types of segregation occurring in the liquid state, and the re- 
sults are found to agree qualitatively with previously published results 
of simulation and of a simple model based on condensation. 

Segregation of hard sphere mixtures has a long history, starting perhaps 
from the classic papers by Wood and Jacobson [1], followed by Lebowitz, 
Rowlinson, and Widom [2,3], extending to the most recent works by vari- 
ous groups [4]. It is still controversial whether or not hard sphere mixtures 
do phase segregate in the absence of gravity [4]. However, in the presence 
of gravity, the situation is rather clear cut. Rosato et al. [5] advanced 
an explanation of the segregation in vertically shaken mixtures of granular 
materials by appealing to geometrical reorganization: during shaking, voids 
opening beneath larger particles are filled more readily by smaller particles. 
This may be contrasted with the apparent "buoyancy" of larger particles [6] , 
or convection driven segregation [7]. But recent molecular dynamics sim- 
ulations of weakly dissipative hard sphere and hard disk granular systems 
under gravity that are in global thermal equilibrium with a heat reservoir 
have shown that a reverse segregation phenomenon occurs as well [7] . These 
results were interpreted in light of a recent proposal that the hard spheres 
undergo condensation transition under gravity [8], which was subsequently 
tested by Molecular Dynamics simulations [10]. 

The density profile in the condensed regime is fairly uniform at the level of 
Enskog approximation [8], or displays oscillatory structure in the weighted 
density functional approximation [9, 11], which is consistent with the ex- 
perimental observation of the formation of crystalline structure with fairly 
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uniform density near the bottom of the shaken granular materials [18]. It 
was argued that the type of segregation (Brazil nut [BN] or reverse Brazil 
nut [RBN]) results from a competition between the system's tendency to re- 
organize geometrically in the ways described by Rosato et al. [5] and the 
tendency of hard spheres under gravity to condense at low enough kinetic 
temperature [8-10]. It was suggested as a qualitative model [7] that in a 
binary hard sphere mixture, the particles belonging to the different species 
would be characterized by different condensation temperatures, so that dur- 
ing quenching from high to low temperature, the species with highest would 
tend to condense at the bottom of the sample first, leading ultimately to 
vertical segregation if this tendency were reahzed. This condensation driven 
tendency toward segregation, however, would be either augmented or op- 
posed by the essentially geometric mechanisms identified by Rosato et al., 
which favor BN segregation only, depending respectively upon whether the 
smaller particle has higher or lower condensation temperature than the larger 
particle. The purpose of this Letter is to present for the first time a theory 
based on the variational principle and directly derive the phase boundary of 
the binary mixtures. 

Theory : We calculate in the local density approximation (LDA). Specif- 
ically, with being the Helmholtz free energy per particle for the ideal gas 
in the absence of gravity, and 11)^x0 being that of the excess component due 
to particle interactions, the binary mixture free energy per area functional is 
given as a function of the densities Pi{v): 



where 1 and 2 are particle indexes, and the total density, p, is the sum of 
the two: p = Pi + P2- This gives the free energy of a columnar sample of 
the system whose transverse cross section has one unit of area. The plane 
z = is the bottom of the container, and for the problem to be strictly 
one dimensional, the size of the container in the transverse direction must 
be infinite. Minimizing F under the global constraints that the number of 
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particles of each species is conserved: 



N,= / drp,(r), t = l,2 (2) 

yields the desired density profiles. 

The forms of ipid are given in [12] and it is standard [12, 13] that ii Z = 
then 

%^ = f (3) 

For the form of Z, we use a generalized Carnahan-Starling equation of state, 
the so called Mansoori, Carnahan, Starling, Leland approximation for a fi- 
nite number of hard sphere species [14]. Defining = f Sr=i P*-^? i^^ = 
0, 3), with Pi = Ni/V, and n being the number of species in the mixture, 
Z is then given by 
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where p = Y17=i Pi- To perform the integration specified in Eq. (|^) for 
n = 2, we recognize that this is an expression for the excess free energy for 
a homogeneous mixture in which pi and p2 are fixed fractions of p, i.e. pi/p 
and P2/P are both constants. Then one readily sees then that the quantities 
may be written as constant multiples of p only, viz. C,a = CaP, which 
defines the constant Ca in terms of pi, P2 and Di and D2. 

Performing the integration specified in Eq. (j^) with the substitution of 
these relations, we find the desired form for ipexc'- 
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(see refs. [15] and [16]). This done, the functional F[pi{z) , p2{z)] defined in 
Eq. ([l|) is then fully specified. Its minimization is accomplished by solving 
for pi{z) and P2{z) the two equations {i = 1,2): 

6F d[p^id\ d^exc . . 

T- = 1 ^ Vexc + P—J h TTLigZ = Xi, (6) 

opi dpi dpi 

where the Lagrange multipliers Aj are introduced to constrain the minimiza- 
tion according to Eq. (g). These equations, though only algebraic, are non- 
linear and highly nontrivial and must therefore be solved numerically via an 
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iterative scheme for Pi{z) at several points z. The two dimensional problem 
is solved analogously. We use the equation of state of hard disks used by 
Jenkins and Mancini [17] 

Z=1 + |1 + 8§1 + ^, (7) 

Ol ^12 ^2 

where cxj = Di/2 and aij = CTj + (Xj, and where 
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Kij = QpiPj^ij9ij- (8) 

The functions Qij are the pair correlations at contact. Following Jenkins and 
Mancini [17], we use 
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where 77, = ^Dfpi, the area fraction of species i, r] = 771 + 772, and R = D1/D2. 
Eq. is integrated, again using the fact that the fractions pi/p are constant 
and then equations analogous to Eq. (^) are then solved iteratively for Pi{z), 
where the derivatives of ipexc are evaluated numerically. 

Results and Discussion: The control parameters that directly enter 
this formulation of the problem are g, T, rrii, Di, and Aj. In practice, for any 
choice of g, T, nii, Di, and pi, the last of which control the average particle 
density, the parameters Aj are tuned so that the integrated density profiles 
yield desired layer numbers pi according to pi = Df~^ f dzpi, where d is the 
dimensionality of the system. To show that the LDA is capable of generat- 
ing reliable results, we present Fig. 0, which shows volume fraction profiles 
r]i{z) = lD^pi{z) for d = ?>, = D2 = 0.001 m, 777i = 1.047 x lO"" kg, 
1712 = 27^1, g = 9.8 m/s^ generated by the LDA (lines) and generated by MD 
simulation (symbols) used in ref. [7]. The layer numbers pi and p2 are both 
nominally 8, but the integrals of the LDA curves disagree slightly with those 
of the MD curves as a result of a systematic binning error in the MD results. 
The temperatures of the two systems differ; Tida/^29D2P2 = 0.0731 and 
Tmd/^29D2P2 = 0.1401. The LDA system is not at all dissipative, so it is 
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Figure 1: A comparison of density profiles calculated from the LDA (lines) 
and from molecular dynamics simulation (symbols). The dashed line and 
squares designate species 1, and the solid line and circles designate species 
2; m2/mi = 2; D2/D1 = 1; Hi ^ ^2 ^ 8. 

reasonable that the different systems give similar profiles only at different 
temperatures, Tmd being necessarily the higher one. We are, however, not 
able to exactly pin point the factor two difference in temperatures. The tem- 
perature difference notwithstanding, the resulting profiles are very similar, 
and indeed exhibit non-BN segregation. Note that in this case, the segrega- 
tion mechanism cannot be geometric, because the particles are of equal size 
but with different mass. 

In all subsequent work reported here, we choose g = 10 m/s^, Di = 
0.001 m, mi = 10^^ kg, and fii = 10 to within less than one tenth percent. 
To investigate the dependence of segregation on mass and diameter ratios, 
we vary D2 and m2 such that D2 > D\ and m2 > mi, find rii{z) at low 
temperature, and quantify the segregation as a function of D2/-D1, m2/mi, 
and T. Our choice to quantify the segregation is simply the ratio of the 
centers of mass of the two species: {zi)/{z2). Specifically, (zi) / {22) < 1 
indicates BN segregation, {zi) / {Z2) = 1 indicates crossover, and (zi) / {Z2) > 1 
indicates RBN segregation. To reduce the number of control parameters and 
to try to make some comparison with the condensation theory [7] which 
originally motivated this work, we choose to explore the lower temperature 
regimes at a fixed reduced temperature. Because the LDA cannot generate 
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information about the solid phase, we are careful to choose temperatures 
at which total volume fractions, ri{z) = rii{z) + r]2{z) have their maximum 
somewhere on the range 0.50 to 0.65 when d = 3, that is, near the density 
of the simple cubic packed single species solid (r/ ^ 0.52), but not as dense 
as the hexagonally packed single species solid {i] ^ 0.74). When d = 2, the 
maximum of ri{z) is kept on the range 0.80 to 0.85 (also between the square 
and triangular packing single species area fractions, approximately 0.79 and 
0.91). Thus, our method explores fluid systems of moderate to high densities 
near the bottom of the sample. We have found that choosing a reduced 
temperature T = T/{m2gD2fJ,2) = 0.0375 {d = 3) and T = T/{m2gD2fi2) = 
0.0300 {d = 2) keeps the maximum of ri{z) on the prescribed ranges over the 
spectrum of diameter and mass ratios we have used. 

Fig. ^ summarizes the results of our calculations for diameter ratios 
D2/D1 = 1.2, 1.5, 2.0, 3.0, 4.0, and a range of mass ratios for both d = 2 
and 3. If from these data we interpolate to find those values of 7712/ mi at 
which crossover occurs {{zi) / {Z2) = 1) for given values of D2/D1, we generate 
the data shown in Fig. ^. Each graph displays three curves. The dashed 
curves with symbols are the crossover curves calculated with the LDA at the 
low temperatures cited above and also at T = 0.1500 for both d = 2 and 
d = 3, so that the temperature dependence may be illustrated. To the left 
of a crossover curve, the segregation is of the BN variety; to the right, it 
is RBN. Thus, the LDA reveals that as reduced temperature increases, an 
increasingly larger region of parameter space gives rise to RBN segregation. 
Indeed, if we assume that in the high temperature limit the density of each 
species is proportional to exp(— m^f/z/T), then for /xi = fi2, one may calculate 
for d = 2 and d = 3 that the ratio of centers of mass is given by 

{zi) m2 



{Z2) mi ' 



(10) 



which indicates the crossover curve is the vertical line — = 1, a result 

mi ' 

consistent with the trend in the LDA results presented in Fig. |^. 

Also appearing on the both graphs in Fig. ^ as solid lines are crossover 
curves predicted from the simple condensation argument. Briefly, these are 
obtained by assuming the ratio of condensation temperatures in the single 
species theory 

Tc2/Tci = m2gD2H2/migDifii 

measures the tendency of species 2 to segregate to the bottom at T < Tc2 
due to the condensation mechanism [7]. 
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Figure 2: Segregation parameter {zi) / (^2) as a function of mass ratio m2/mi 
for hard spheres (left) and hard disks (right). In both graphs Hi = IJ.2 = 10; 
D2/D1 = 1.2 (circle), 1.5 (square), 2.0 (diamond), 3.0 (up triangle), 4.0 
(down triangle), fg^ = 0.0375; fad = 0.0300. 
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Figure 3: Curves through parameter space along which the segregation pa- 
rameter (zi) I (^2) = 1; indicating the crossover from BN (left) to RBN (right) 
segregation, for (i = 3 (upper) and for d = 2. The dashed curves marked 
by circles are from linear interpolation of the data shown in Fig. ^ and are 
for the low temperatures given in the text. The dashed curves marked by 
squares are for higher reduced temperature (T = 0.1500 in both cases). 
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If we choose to equate this with Rosato et aL's [5] control parameter for 
the reorganization, {D2/D1Y, we can obtain the crossover curve (The sohd 
hues in Fig. 3). We see fair agreement with our theoretical prediction at 
low reduced temperature for (i = 3, but significant disagreement appears at 
1712/1711 ^ 7; the crossover curve seems practically linear, in disagreement 
with the scaling suggested by the condensation argument. For d = 2, the 
low temperature LDA boundary is linear, in agreement with the scaling sug- 
gested by the condensation argument, but that argument of course cannot 
predict the slope. Moreover, the LDA computed boundaries are necessarily 
temperature dependent, but the condensation theory [7] does not predict a 
temperature dependence. Thus, although the condensation theory has been 
useful in directing us to look for certain dependencies and, indeed, for the 
phenomenon of RBN segregation itself, its utility, except as a heuristic, is 
somewhat circumscribed. However, the equilibrium approach of minimizing 
Helmholtz free energy seems a promising theoretical alternative, as it is ca- 
pable of generating solutions in good agreement with those obtained from 
MD simulation, as in Fig. |I|. 
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